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I. INTRODUCTION 

Since the mechanical exfoliation of graphene^, research 
in understanding the properties of two-dimensional (2D) 
crystals has increased in many folds. Atomically thin 
single-layered materials obtained from transition metal 
dichalcogenides^, boron nitride^, bismuth^, etc. are be¬ 
ing extensively studied for applications as electronic and 
photoelectronic devices. Few-layered black phosphorous 
(BP) is a recent addition to the list of graphene-inspired 
materials^”^. Apart from having a sizeable band gap 
which can be tuned by the manipulation of the number of 
layers, the atomic structure of BP is highly anisotropic 
which leads to high asymmetry of the electronic band 
structure even for few layers. In particular, a single layer 
of BP or phosphorene is attracting most of the attention. 
The peculiar anisotropic nature of the band gap distin¬ 
guishes phosphorene from other 2D crystals, increasing 
its potential functionality. 

Excitons are a bound state of an electron and an hole 
and play an important role in the optical properties of 
the material. Understanding the nature of excitons and 
their dependence on the electronic structure of the host 
material is critical and lends a deeper perspective into 
the many-body physics involved in 2D crystals. The 2D 
nature of the polarizability of these crystals introduces 
an important length scale (screening length) vq. For dis¬ 
tances between charges in the crystal plane, r, greater 
than ro the electron-hole binding potential behaves like 
in a 3D system i.e., it goes as ^ 1/r. However, for the 
case where r is less than ro the potential is 2D-like, i.e., 
logarithmic. This behavior makes excitons in 2D crystals 
different from their 3D counterparts^’ 

Since most common 2D crystals are isotropic, the effect 
of anisotropy on the optical properties of these materials 
has remained essentially unexplored. The appearance of 
phosphorene has, however, changed this view and recent 
works address this issue from an analytical^, numerical^ 


and first-principles^ standpoints. Here we give a detailed 
account of a variational approach, introduced by us in 
Ref. 8, to the calculation of the exciton binding energy 
Ex in anisotropic 2D crystals. Several analytical expres¬ 
sions are derived in certain limits of the 2D interaction 
potential. The accuracy of our analytical expressions for 
Ex is tested against both variational and exact numeri¬ 
cal solutions to the actual 2D potential, finding excellent 
agreement in a wide and experimentally relevant range 
of screening lengths. In particular, the value of Ex for 
phosphorene, as obtained from our analytical expression, 
compares almost exactly to the numerical results. Fur¬ 
thermore, this value nicely agrees with the recently re¬ 
ported experimental result^^. 


The present work is divided as follows. In Sec. H we 
review the form and limiting behaviour of the Coulomb 
interaction potential for charged particles in 2D systems. 
In Sec. HI we present our variational approach based on 
an anisotropic exciton wavefunction. We first present the 
analytical result for Ex in the limiting case where the 2D 
potential reduces to the standard 3D Coulomb potential 
^ 1/r for isotropic 2D systems to later introduce the 
anisotropy and re-derive the binding energy for this case. 
In the same manner we derive analytical expressions for 
the isotropic and anisotropic binding energies in the op¬ 
posite limit where the 2D interaction potential behaves 
logarithmically. We also compare our analytical expres¬ 
sions with the numerically solved variational problem as 
well as with the exact numerical solution. In Sec. IV 
we propose an alternative variational approach based on 
gaussian orbitals. In Sec. V, after computing the 2D 
polarizability with density functional theory (DFT), our 
analytical approach is applied to the case phosphorene. 
Finally we present our conclusions in Sec. VI. 
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II. BINDING PARTICLE-HOLE POTENTIALS 
IN 2D 


As originally derived by Keldysh^^, the Coulomb po¬ 
tential energy created by a point charge at the origin that 
electrons feel in 2D layers follows the expression: 


h^ 2 D(r) 


Seoero 




( 1 ) 


where ro = de/(ei + 62 ) and e = (ei + e 2 )/ 2 . Here d is 
the thickness of the 2D material, e is its bulk dielectric 
constant, and ei and 62 are the dielectric constants of 
the surrounding media, typically substrate and vacuum. 
Here ro plays the role of a screening length and sets the 
boundary between two different behaviours of the poten¬ 
tial. For r < ro the potential diverges logarithmically, 
as if created by line charges. In this limit, the potential 
takes the simplified form also given by Keldysh^^: 


V 2 D(r < ro) 


1 

dTreoe ro 



( 2 ) 


where 7 is the Euler constant. For r > tq the poten¬ 
tial becomes the standard Coulomb potential created by 
point charges which decays as 1 /r: 


V 2 B{r > ro) 


1 

47 reoe r 


(3) 


A very good approximation to the Keldysh potential, 
fairly accurate in both limits and simpler to use, was 
introduced by Cudazzo et al.^\ 


>” (;rT7^) + ■ 

( 4 ) 

It is interesting to compare these four expressions as 
a function of the distance r in a range of several orders 
of magnitude both above and below ro- We present such 
a comparison in Fig. 1. There it can seen the range of 
validity of each approximation, the Cudazzo et al. ex¬ 
pression being remarkably accurate for all distances. 


^ 2 d(^) 


Aneoe vq 


III. VARIATIONAL WAVEFUNCTION 
APPROACH 


For generic 2D crystals with electrons and holes pre¬ 
senting anisotropic effective masses, we 

consider variational solutions for the exciton wave func¬ 
tion of the type^^: 



where A is the variational anisotropy scaling factor relat¬ 
ing the exciton extension along the x direction, (which 
is also a variational parameter) and the one along the y 



FIG. 1: (color online). Keldysh, 3D Coulomb, 
logarithmic, and Cudazzo et al. potentials in log-log 
scale in a range of distances spanning several decades 
around Vq. 


direction {ay = Aa^,). With this variational wavefunc- 
tion we can evaluate the expectation value of the kinetic 
energy: 


-F'kin {(^x 5 ^) 
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where /j,x and jXy are the reduced effective masses, 
m®m^/(m® +m^), along x and y directions, respectively. 
The expectation value of the potential energy is given by 



V 2 B{x,y)(f>{x,yfdxdy 


(6) 


and the variational exciton binding energy is obtained 
from the addition of these two quantities. 


-^x(^a:5 4^) -^kin “1“ -^pot • (7) 

Upon minimization with respect to and A, one obtains 
the optimal parameters defining the extension and shape 
of the exciton and the actual binding energy Ex- Results 
from three minimization procedures, one analytical and 
two numerical, are presented in next section. 


A. Analytical Results 

The integral for the potential energy in Eq. ( 6 ) turns 
out to be too difficult for an exact variational analytical 
solution. The main goal of this section is to make use 
of the asymptotic behaviour of the Keldysh potential to 
get analytical expressions for Ex in the limits r ^ ro 
and r <C ro, namely, valid for large and small excitons, 
respectively. 
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1. r ^ tq limit 


We begin by evaluating Ex in the isotropic case (A = 1 , 
= a), considering only the long-range behaviour 
of the Keldysh potential (see Eq. 3). The contribution 
of the potential energy to Ex is given in this limit by 


E, 


pot 
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(8) 


Now minimizing Ex (cl) with respect to the variational 
exciton radius one obtains 


Ex = — 


1 

47reo6 a ’ 


(9) 


where the minimal exciton radius a is given by 


a = 


a^em 

2 /i ’ 


( 10 ) 


and m and an = are the free electron mass and 

the Bohr radius, respectively. 

For the anisotropic case (A 7 ^ 1) the exciton extension 
along the x-direct ion is now given by 


We now define 


Ie{\)= {>? - l)dl{\)/d\ + A/(A) 


(17) 


where E is the complete elliptic integral of second kind. 
We can see that the minimal A, A, satisfies in general the 
following equation: 


Ai. 3 /s(A)-A7(A) 

which has no analytical solution for A. However, it can 
be shown^^ that for A < 1 


/ \ 1/3 

Finally, notice that the results obtained in this subsection 
will be valid as long as the exciton extension in both x 
and y directions is much larger than rp. The consistency 
of this approximation for given experimental parameters 
(ro, My, and e) has to be checked a posteriori. 


2. r ^ ro limit 


- aoem f 1 1 ^ 1 

dxW = -^ I-h 
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In the previous expression we find a function of A defined 
through the elliptic integral 


1 


■ 


A 


+ (A^ — 1 ) cos^ 0 

( 12 ) 
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where the function K is the complete elliptic integral of 
the first kind. Defining now ^xy as 
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l^x X^ fJj 
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/(A), (13) 


the exciton extension along the x axis can now be written 
as 


.(A) = 


ao em 


2 /i'xy(A) 

The exciton extension along the y direction is thus 


(14) 


. ... aoem / I 1 \ A ao emX 
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and the A-dependent binding energy of the exciton now 
becomes 


Ex{X) = - 


AweoedxiX)' 


(16) 


As r ^ 0 the logarithmic behaviour of the Keldysh 
potential dominates. The potential energy in 2D takes 
now the form given in Eq. (2). In the isotropic case the 
exciton radius is now given by 


a = 


em 




-doTo 


and the binding energy of the exciton is 


Ex = 


1 


Attcoc ro 


In 
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( 20 ) 


( 21 ) 


For an anisotropic system the A-dependent exciton ex¬ 
tension along the x direction is given by 


- /AN / em f I 1 

«.(A) = (^- + ^ 

Using now a different definition for jixy 


( 22 ) 
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the exciton x-extension now becomes 


= \ dor0 


(24) 

y Mx7/(A) 

Again, taking into account that ay = Aa^,, the A- 
dependent minimal exciton extension along the y direc¬ 
tion is 


j{X) = Waoro 


emX^ 

Mxy (A) 


(25) 



Note th.8jt CLx^^'X’) l^y) -) l^x) • 

Finally we obtain the exciton energy for this case: 
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where the minimal A is 

for all jiix and jiy. Again, notice that this result will be 
valid as long as the x and y minimal extensions of the 
excitonic wave function are small compared to tq. 

Note that Eq. (26) can be written in a more symmet¬ 
rical way as a function of both and d^, 


Ex{\) 


1 3 ^ /da^(A) + dy(A) 

dTreoero ^2 V Svq 


and that it is also symmetrical under exchange of jXx and 
as it should be. However, we find that the binding 
energy is not only a function of yxIfJ^y^ but depends on 
both their values. 

Finally, for completeness, we present an analytical ex¬ 
pression for the exciton binding energy using the Cudazzo 
potential in the isotropic case: 


Ex 


do 


+ 4(7-ln(2))- 


ro 
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dro 
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a 
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(28) 


where Ei is the exponential integral function. We have 
only been able to obtain a working analytical expression 
for d (too cumbersome to be shown here) in the limit 
^0 ^ ^^0 where the above expression is actually useful. 


B. Numerical Optimization and exact solution 

To validate and test the accuracy of the limiting ana¬ 
lytical expressions given in the previous section, we now 
use the wavefunction in Eq. (5) to numerically com¬ 
pute the potential energy given by Eq. (6) for the exact 
Keldysh potential. We also solve, numerically as well, the 
2D Schrodinger equation for the same potential, which 
will give us the exact value of Ex (down to the required 
numerical precision). Exciton binding energies for the 
isotropic case (A = 1) are presented in Fig. 2 as a func¬ 
tion of the screening length vq. For comparison’s sake, 
we take e = 1, i.e., the 2D crystal is suspended in vac¬ 
uum, and yx = yy = m/2. Thus, according to Eq. (9), 
Ex = —2 Ryd (where Ryd is the Rydberg energy 13.6 
eV) for ro = 0. The numerical variational result com¬ 
pares very well with the exact numerical value in the large 
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FIG. 2: (color online). Binding energies (in Ryd) for 
the isotropic exciton computed in different ways: 
Numerical optimization of the variational wavefunction 
in Eq. 5 (solid line), numerical solution of the 
Schrodinger equation (circles), analytical variational 
expression valid for large values of ro (Eq. 21) (dashed 
line), and expression in Eq. 28 obtained for the 
Cudazzo et al. potential (dotted-dashed line) 


range of explored screening lengths. For ro ^ no the an¬ 
alytical solution in Eq. (21) works fairly well. There, 
the size of the exciton is smaller than ro and the 1/r 
contribution to the Keldysh potential is negligible. As 
expected, the analytical solution starts to fail as ro ^ no 
since there the size of the exciton becomes comparable 
to ro and the long-range 1/r contribution to the Keldysh 
potential becomes dominant. (One should keep in mind 
that the limit of validity of the analytical result, as shown 
in Eq. (20), depends on the values of e and y.) We also 
compare with the result given by Eq. (28), obtained us¬ 
ing the approximate expression to the potential in Eq. 
(4). This expression, although not as friendly as the pre¬ 
vious one, extends the limit of validity of our analytical 
results down to ro ~ ao- 


The results for the anisotropic case are presented in 
Fig. 3 for ro = 20ao and ro = 400ao as a function of the 
anisotropy ratio ^ with yx = m/2 (notice a difference of 
one order of magnitude in the energy scales of each plot). 
Note that these curves would be identical if plotted as a 
function of — with Uy = m/2. Once again there is close 

fly y 

agreement between the analytical solution [Eq. (26)], the 
numerical optimization, and the exact numerical solution 
for large ro, while for the smaller value, the analytical 
solution visibly deviates from the other two. 


An important prediction of our analytical results is the 
relation between the anisotropy in the exciton extension 

and the effective masses: A = ^ {^) ’ which 


becomes exact in the limit of small excitons. To test this 
relation we fitted the optimal value of the variational 
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as obtained from the analytical and the numerical 
approaches for vq = 20ao (a) and vq = 400ao (b). 


parameter A to the law 


A = C(ro)(^ 


a(ro) 


(29) 


FIG. 5: (Color online). Exact (dashed lines) and 
variational (solid line) wave functions for 
ro = 0.1,1,10,100. The x-axis is rescaled with the 
variational extension of the exciton and the amplitude 

with the normalization constant of the variational wave 

/ \l/2 

function A = 


for a large range of vq. The results of this fit are presented 
in Fig. 4. They confirm our analytical results and recover 
the exact 1/3 exponent in the limit of large vq. 

We finally provide a comparison between the exact 
and variational wave functions for several values of ro 
in the isotropic limit (see Fig. 5). Note that the distance 
is rescaled with the optimal radius and the amplitude 
of the wave function with the normalization constant 
( 2 

A = ( 1 .This representation illustrates to what 

extent the exact and variational wave functions satisfy 


similar scaling relations. Note that at r = 0 the exact 
wave functions do not show the prominent cusp of a Is 
Slater-type orbital. This softened behavior at the origin 
suggests than a combination of gaussian functions may 
capture more accurately this feature of the wave function, 
as shown in next section. 
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IV. GAUSSIAN-BASIS VARIATIONAL 
METHOD 

We have found that in the limit of very small ro, the 
binding energy is very sensitive to small changes in ro- 
Furthermore, the numerical solution of the Schrodinger 
equation requires a very fine mesh to reproduce the 
bound state in such a limit. On the other hand, the 
analytical result found in the 1/r limit of the potential 
constitutes an isolated point an thus cannot be easily 
extended to small but finite values of ro- It is there¬ 
fore interesting to find an alternative numerical method 
to study anisotropic excitons in the limit of small ro- 
Moreover, as we have presented in Fig. 5, the behavior 
of the exciton exact wave functions for different values 
of the screening length ro resembles more a Is Gaussian 
than a Slater-type orbital. Gaussian-type orbitals (GTO) 
are very efficient basis sets used intensively in quantum 
chemistry and solid state calculations. 

The gaussian basis functions {Xp} follow the expres¬ 
sion: 



FIG. 6: (Color online). Ground state energies obtained 
using the gaussian variational method (continuous lines) 
and the numerical solution of Schrodinger equation for 
ro = lOao (circles) as a function of the effective mass 
ratio i^y/i^x for = m. 


= (30) 

The index p is an integer that here has been chosen to 
run from 1 to 4 and the exponents a are coefficients that 
have to be optimized to minimize the ground state energy 
obtained by the variational method. Unlike the conven¬ 
tional gaussian approach, the anisotropy of the problem 
introduces two coefficients a'^ per GTO. We limit the 
variational freedom assuming that the anisotropy is iden¬ 
tical for the four basis functions: 


al = KaP. (31) 

In this equation, is a constant that does not depend 
on p so that we reduce the number of exponents that we 
have to optimize from eight to four. The variational wave 
function in the GTO basis is given by 

4 

0g(^5^) = ^ ^ ^pXp' (32) 


For fixed values of the energy is computed by gen¬ 
eralized diagnonalization of a 4x4 matrix. The matrix 
elements of the kinetic energy are 


^kin ^ ^ M « , 1 \ 

Mpq + otl pyal + alJ' 

The matrix elements of the potential energy are 
puted by numerical integration, 

^pq — // ^2D (^5 y)XpXqdxdy ^ 

and the overlap matrix is 


(33) 

com- 


(34) 


( 35 ) 


where Mpq = y/-h al). 

The binding energy and the optimal wave function are 
obtained by minimizing numerically the energy with re¬ 
spect to the five variational parameters. Efficient opti¬ 
mization of the energy requires a careful choice of the ini¬ 
tial guess for the values of the exponents . In our case, 
we choose the optimal values for a Is orbital of a “2D hy¬ 
drogen atom”, taking = py = m and ro approaching 
zero. Once these optimal exponents are obtained, ro is 
changed slightly and the problem is solved again, using 
this time the exponents a obtained in the previous step. 
The procedure continues until the ground state energy 
for the desired ro is reached. 

For example, fixing ro = lOao and px = the effec¬ 
tive mass along the y axis py is varied from 1 to 40m. 
In this case the initial guess for the exponents is the last 
set of coefficients a obtained when changing ro- The 
result for the ground state energy of this calculation is 
shown in Fig. 6 in comparison to the numerical solution 
of the Schrodinger equation with the Keldysh potential. 
They match perfectly. Two other values of ro/ao ob¬ 
tained with the same gaussian-basis variational method 
are also shown. 

In Fig. 7, the length scales = \/\f^ of the whole 
set of optimal GTO’s for ro = lOao are plotted vs PyjPx 
for px = m. K is also plotted in Fig. 8 against the same 
quantity. In Fig. 9 we show a log-log representation of 
hi versus the asymmetry ratio for different values of ro 
where a linear fitting has been made. 


V. APPLICATION TO PHOSPHORENE 

A paramount example of an anisotropic 2D crystal is 
phosphorene^’^’^, where the effective masses along x and 
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y directions can differ by even an order of magnitude. 
From the start we have chosen to express the 2D potential 
constant in terms of the bulk dielectric constant e and 
the effective thickness of the 2D crystal d (see Eq. 1). 
Equivalent expressions for the 2D potential, which rely 
on the evaluation of the actual 2D polarizability of the 
2D crystal, y, have recently been proposed^ These are 
probably more appropriate for actual crystals, although 
it has also been shown that Eq. (1) works well as long as 
e is taken as the in-plane component of the bulk dielectric 
tensor^^ of the 3D crystal. Here we will compare both 
possibilities. 

As shown in Ref. 9, the screening length vq depends 
on the polarizability y as ro = 27ry. The polarizability 
for 2D materials can be computed using the expression 


EIG. 7: (Color online). Coefficients = 1/v^^ for 
ro = lOao obtained versus [iyjyx foi* = m. 


e(L) — 1 + 


47ry 

L 


(36) 



EIC. 8: (Color online). Coefficient n = a^ja^ for 
ro = lOao obtained versus [iyjyx for Mcc = m. 



EIC. 9: (Color online). Log-log representation of the 
coefficient k — obtained for different values of ro 

(circles) and its linear fitting (continuous line) versus 
jiy/yx for = rn. 


where L is the distance between layers in a 3D lay¬ 
ered structure. As can be seen, the dielectric func¬ 
tion e tends to unity as the inter layer distance L 
tends to infinity. We have computed the dielectric 
function at different inter layer distances within the 
density functional theory framework using the Perdew- 
Burke-Ernzerhof (PBE) functional^^ and norm conserv¬ 
ing Troullier-Martins (TM) pseudopotentials as available 
in the SIESTA package^^. The atomic and electronic 
structures have been duly converged on all parameters. 
SIESTA calculates the imaginary part of the dielectric 
function from which the real part of it is obtained using 
the Kramers-Kronig relations. In order to account for the 
under-estimated band gap, the scissor approximation, as 
implemented in SIESTA, has been utilized. The scissor 
shift of 1.2505 eV was made to match our previously re¬ 
ported band gap value of 2.15 eV^. While more elegant 
approaches to the gap problem of phosphorene have been 
reported in the literature^, the scissor approximation suf¬ 
fices to our purpose here. 

Using the plane-averaged static dielectric function cal¬ 
culated with SIESTA (see Eig. 10), a value in the vicinity 
of y of 3.8 A is obtained. This value for y yields a screen¬ 
ing length of ro = 23.2 A. Erom the numerical variational 
solution we obtain an exciton extension of d^ = 11.7 A 
and dy = b.9 A along the x and y directions, respectively. 
Since these values are smaller than ro, the use of the an¬ 
alytical expressions obtained in the logarithmic limit of 
the potential is justified. Also this was expected from 
Eig. 2 and, in particular, from the comparison shown in 
Eig. 3 for anisotropic cases. There it can be seen that al¬ 
ready for ro = 20ao the deviation between the analytical 
result and the numerical ones is less than 10% for a ratio 
~ 7 (which corresponds to phosphorene). Using 
now Eqs. (22)-(27), we obtain an exciton binding energy 
for phosphorene in vacuum of Ex = 0.61 eV, while the 
numerical variational value is Ex = 0.78 eV. This result 
is remarkably close to a recently reported experimental 
value of ~ 0.9 eV. The agreement is somewhat surprising 
since this has been measured for phosphorene on a SiO 
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FIG. 10: (Color online), (a) Real part of the x and y 
components of the dielectric constant (evaluated at zero 
frequency) for different values of the interlayer distance 
L in a 3D black phosphorous structure, (b) 
Polarizability as obtained from Eq. 36 for different 
values of L after averaging on the plane. 




substrate^^. A more recent experiment, however, reports 
a smaller value for Ex^ which is maybe more expected 
due to the screening of the substrate^^. 


VI. CONCLUSIONS 

We have shown that a variational approach to the exci- 
ton binding energy in anisotropic 2D crystals can give ex¬ 
cellent results when compared to numerical approaches. 
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FIG. 11: (Color online). Exciton binding energy (in eV) 
as a function of the thickness when all the other 
parameters (effective mass, bulk dielectric constant, and 
screening length) to correspond to phosphorene. The 
horizontal line marks the value obtained using the 
actual 2D polarizability of phosphorene. 


Furthermore, we have studied the range of validity of an¬ 
alytical solutions to the variational approach and found 
that these can give highly satisfactory results in a range 
of values of screening lengths which is relevant for actual 
2D crystals such as phosphorene. We have computed 
the exciton binding energy in this case and found a very 
good agreement with a recently reported experimental 
result^^. As long as the screening length to exciton size 
is large, our analytical results can be trivially used to 
predict the exciton binding energy of any 2D crystal. 
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